{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using SAGE\n",
    "\n",
    "This notebook shows how to use SAGE to identify author-specific keywords, using a dataset in the format of word count statistics. The dataset is from [Ted Underwood](https://tedunderwood.com/open-data/), and the relevant subsets are included as [releases](https://github.com/jacobeisenstein/SAGE/releases/tag/data) in this repository.\n",
    "\n",
    "If you use SAGE, please cite [the paper](http://www.icml-2011.org/papers/534_icmlpaper.pdf):\n",
    "\n",
    "```Eisenstein, Jacob, Amr Ahmed, and Eric P. Xing. \"Sparse Additive Generative Models of Text.\" Proceedings of the 28th International Conference on Machine Learning (ICML-11). 2011.```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import sage\n",
    "from collections import Counter\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we build word count dictionaries from each file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def getCountDict(filename):\n",
    "    with open(filename) as fin:\n",
    "        return {word:int(count) for word,count in [line.rstrip().split() for line in fin.readlines()]}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# counts for author Lydia Maria Child\n",
    "child_counts = getCountDict('underwood-child-counts.tsv')\n",
    "# counts for all 1840s letters in the corpus \n",
    "base_counts = getCountDict('underwood-1840s-let-counts.tsv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Build a vocabulary of the most common terms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "vocab = [word for word,count in Counter(child_counts).most_common(5000)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convert the counts into [numpy](http://www.numpy.org/) arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "x_child = np.array([child_counts[word] for word in vocab])\n",
    "x_base = np.array([base_counts[word] for word in vocab]) + 1."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the base log-probabilities of each word"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "mu = np.log(x_base) - np.log(x_base.sum())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run SAGE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "eta = sage.estimate(x_child,mu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print words that are especially frequent in the writing of Lydia Maria Child, compared to the baseline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['soul', 'sea', 'heart', 'beauty', 'outward', 'graceful', 'mother', 'souls', 'forth', 'sweet']\n"
     ]
    }
   ],
   "source": [
    "print sage.topK(eta,vocab)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print words that are especially infrequent in her writing, compared to the baseline."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['america', 'xx', 'extent', 'democratic', 'states', 'population', 'parties', 'falls', 'americans', 'timber']\n"
     ]
    }
   ],
   "source": [
    "print sage.topK(-eta,vocab)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A histograph of SAGE coefficients. The spike at zero is because SAGE is biased towards sparsity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAEACAYAAABfxaZOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFh5JREFUeJzt3X+Q5HWd3/Hna9lAOD03qHe71q6CCt6hKQOmBFKk6hoN\nuEAVmKsc2btL8UNNaYl3VOVyYUHjDmrVCXV44F1Z/CFaYIVbicYDDcJqoE2ZHAu6rIK7wJoEXFbY\nu1KXO/CKQnjnj/4MtsOM0/Njp2fn+3xUTe23P/35dr97dubVn35/vz2dqkKS1A2rxl2AJGnpGPqS\n1CGGviR1iKEvSR1i6EtShxj6ktQhI4d+klVJdiS5tV0+JsndSR5O8pdJVrfxw5NsTbInyV8nec3Q\nbVzWxncnOWPxH44k6ZeZy0r/EmDX0OUrgaur6g3AAeDdbfzdwI+r6jjgGuAqgCRvBM4DjgfOBD6V\nJAsrX5I0FyOFfpINwFnAp4eG3wZ8sW3fALyzbZ/bLgN8oc0DOAfYWlU/q6pHgD3ASfOuXJI0Z6Ou\n9P8M+GOgAJK8AvhJVT3frn8MWN+21wN7AarqOeDJJC8fHm/2De0jSVoCs4Z+krOB/VW1Exhux4za\nmrGFI0nLxOoR5pwKnJPkLOBI4FeBa4E1SVa11f4GBit32r+vBn6Y5DDgZVX14yST45OG93lBEv8Y\nkCTNQ1XNusiedaVfVZdX1Wuq6nXAJuDOqvp3wF3A77RpFwC3tO1b22Xa9XcOjW9qZ/e8FjgWuGeG\n+1x2X1u2bBl7DdZkTV2sy5pG+xrVKCv9mWwGtib5KHAfcH0bvx74XJI9wI8YPFFQVbuS3MzgDKBn\ngffXXCqVJC3YnEK/qr4BfKNt/z/g5GnmPMPg1Mzp9v8T4E/mXqYkaTH4jtwR9Xq9cZfwItY0Gmsa\n3XKsy5oWV5ZbhyWJXR9JmqMk1GIcyJUkrRyGviR1iKEvSR1i6EtShxj6ktQhhr4kdYihL0kdYuhL\nUocY+pLUIYa+JHWIoS9JHWLoS1KHGPqS1CGGviR1iKGvQ8a6dceQZN5f69YdM+6HII2df09fh4wk\nwEJ+NjKnzxKVDiX+PX1J0ovMGvpJjkiyPcl9Se5PsqWNfzbJ/23jO5K8eWifTybZk2RnkhOGxi9I\n8nCSh5Kcf3AekiRpJrN+MHpVPZPktKr6aZLDgP+V5PZ29X+sqv82PD/JmcDrq+q4JCcD1wGnJDkK\n+DDwFiDAt5PcUlVPLuojkiTNaKT2TlX9tG0eweCJ4vl2ebr+0bnAjW2/7cCaJGuBdwDbqurJqjoA\nbAM2LqB2SdIcjRT6SVYluQ94AvhaVd3brvpYa+FcneQftbH1wN6h3R9rY1PH97UxSdISGXWl/3xV\nnQhsAE5K8kZgc1UdD7wVeAVw6Qy7z3o0WZK0NGbt6Q+rqr9L0gc2VtUn2tizST4L/FGbtg949dBu\nG9rYPqA3Zfyu6e5nYmLihe1er0ev15tumiR1Vr/fp9/vz3m/Wc/TT/JK4NmqejLJkcAdwMeBHVX1\nRAYnT38C+IequjzJWcDFVXV2klOAa6pq8kDutxgcyF3Vtv956+8P35/n6WtanqcvzWzU8/RHWem/\nCrghySoGYf35qrotyf9oTwgBdgLvA2jXnZXk+8DTwEVt/CdJPsog7Au4YmrgS5IOLt+Rq0OGK31p\nZr4jV5L0Ioa+JHWIoS9JHWLoS1KHGPqS1CGGviR1iKEvSR1i6EtShxj6ktQhhr4kdYihL0kdYuhL\nUocY+pLUIYa+JHWIoS9JHWLoS1KHGPqS1CGGviR1iKEvSR0ya+gnOSLJ9iT3Jbk/yZY2fkySu5M8\nnOQvk6xu44cn2ZpkT5K/TvKaodu6rI3vTnLGwXtYkqTpzBr6VfUMcFpVnQicAJyZ5GTgSuDqqnoD\ncAB4d9vl3cCPq+o44BrgKoAkbwTOA44HzgQ+lcEnXUuSlshI7Z2q+mnbPAJYDRRwGvDFNn4D8M62\nfW67DPAF4G1t+xxga1X9rKoeAfYAJy2keEnS3IwU+klWJbkPeAL4GvB/gANV9Xyb8hiwvm2vB/YC\nVNVzwJNJXj483uwb2keStARWjzKphfuJSV4GfAn4zTncx5xbOBMTEy9s93o9er3eXG9Ckla0fr9P\nv9+f836pqrntkPxn4B+A/wSsq6rnk5wCbKmqM5Pc3ra3JzkMeLyqfj3JZqCq6sp2Oy/Mm3L7Ndea\n1A2DQ0AL+dkI/mxppUpCVc26yB7l7J1XJlnTto8ETgd2AXcBv9OmXQDc0rZvbZdp1985NL6pnd3z\nWuBY4J7RHo4kaTGM0t55FXBDklUMniQ+X1W3JdkNbE3yUeA+4Po2/3rgc0n2AD8CNgFU1a4kNzN4\nwngWeL9LeklaWnNu7xxstnc0E9s70swWrb0jSVo5DH1J6hBDX5I6xNCXpA4x9CWpQwx9SeoQQ1+S\nOsTQl6QOMfQlqUMMfUnqEENfkjrE0JekDjH0JalDDH1J6hBDX5I6xNCXpA4x9CWpQwx9SeoQQ1+S\nOmTW0E+yIcmdSb6X5P4kf9DGtyR5LMmO9rVxaJ/LkuxJsjvJGUPjG5M8mOThJJcenIckSZrJrB+M\nnmQdsK6qdiZ5KfBt4Fzg3wJ/X1WfmDL/eOAm4K3ABuDrwHFAgIeBtwM/BO4FNlXVg1P294PRNS0/\nGF2a2agfjL56tglV9QTwRNt+KsluYP3k/Uyzy7nA1qr6GfBIkj3ASW3unqp6tBW4tc19cJrbkCQd\nBHPq6Sc5BjgB2N6GLk6yM8mnk6xpY+uBvUO77WtjU8cf4+dPHpKkJTDrSn9Sa+18Abikrfg/BXyk\nqirJx4CrgfcsRlETExMvbPd6PXq93mLcrCStGP1+n36/P+f9Zu3pAyRZDXwF+GpVXTvN9UcDX66q\nNyfZDFRVXdmuux3YwqC9M1FVG9v4L8wbui17+pqWPX1pZqP29Edt73wG2DUc+O0A76TfBh5o27cC\nm5IcnuS1wLHAPQwO3B6b5OgkhwOb2lxJ0hKZtb2T5FTg94H7k9zHYKl1OfB7SU4AngceAd4LUFW7\nktwM7AKeBd7flu7PJfkAsI3Bk831VbV78R+SJGkmI7V3lpLtHc3E9o40s8Vu70iSVgBDX5I6xNCX\npA4x9CWpQwx9SeoQQ1+SOsTQl6QOMfQlqUMMfUnqEENfkjrE0JekDjH0JalDDH1J6hBDX5I6xNCX\npA4x9CWpQwx9SeoQQ1+SOsTQl6QOmTX0k2xIcmeS7yW5P8kftvGjkmxL8lCSO5KsGdrnk0n2JNnZ\nPjx9cvyCJA+3fc4/OA9JkjSTWT8YPck6YF1V7UzyUuDbwLnARcCPquqqJJcCR1XV5iRnAh+oqrOT\nnAxcW1WnJDkK+BbwFiDtdt5SVU9OuT8/GF3T8oPRpZkt2gejV9UTVbWzbT8F7AY2MAj+G9q0G9pl\n2r83tvnbgTVJ1gLvALZV1ZNVdQDYBmyc06OSJC3InHr6SY4BTgDuBtZW1X4YPDEAa9u09cDeod0e\na2NTx/e1MUnSElk96sTW2vkCcElVPZVk6uvkmV43z/pyY6qJiYkXtnu9Hr1eb643IUkrWr/fp9/v\nz3m/WXv6AElWA18BvlpV17ax3UCvqva3vv9dVXV8kuva9ufbvAeB3wJOa/Pf18Z/Yd7QfdnT17Ts\n6UszW7SefvMZYNdk4De3Ahe27QuBW4bGz29FnAIcaG2gO4DTk6xpB3VPb2OSpCUyytk7pwL/E7if\nwTKrgMuBe4CbgVcDjwLntQO0JPkLBgdpnwYuqqodbfxC4IPtNj5WVTdOc3+u9DUtV/rSzEZd6Y/U\n3llKhr5mYuhLM1vs9o4kaQUw9CWpQwx9SeoQQ1+SOsTQl6QOMfQlqUMMfUnqEENfkjrE0JekDjH0\nJalDDH1J6hBDX5I6xNCXpA4x9CWpQwx9SeoQQ1+SOsTQl6QOMfQlqUMMfUnqkFlDP8n1SfYn+e7Q\n2JYkjyXZ0b42Dl13WZI9SXYnOWNofGOSB5M8nOTSxX8okqTZzPrB6En+JfAUcGNVvbmNbQH+vqo+\nMWXu8cBNwFuBDcDXgeOAAA8Dbwd+CNwLbKqqB6e5Pz8YXdPyg9GlmY36weirZ5tQVd9McvR09zHN\n2LnA1qr6GfBIkj3ASW3unqp6tBW3tc19UehLkg6ehfT0L06yM8mnk6xpY+uBvUNz9rWxqeOPtTFJ\n0hKadaU/g08BH6mqSvIx4GrgPYtV1MTExAvbvV6PXq+3WDctSStCv9+n3+/Peb9Ze/oArb3z5cme\n/kzXJdkMVFVd2a67HdjCoL0zUVUb2/gvzJtye/b0NS17+tLMRu3pj9reCUM9/CTrhq77beCBtn0r\nsCnJ4UleCxwL3MPgwO2xSY5Ocjiwqc2VJC2hWds7SW4CesArkvyAwcr9tCQnAM8DjwDvBaiqXUlu\nBnYBzwLvb8v255J8ANjG4Inm+qravfgPR5L0y4zU3llKtnc0E9s70swWu70jSVoBDH1J6hBDX5I6\nxNCXpA4x9CWpQwx9SeoQQ1+SOsTQl6QOMfQlqUMMfUnqEENfkjrE0JekDjH0JalDDH1J6hBDX5I6\nxNCXpA4x9CWpQwx9SeoQQ1+SOmTW0E9yfZL9Sb47NHZUkm1JHkpyR5I1Q9d9MsmeJDvbh6dPjl+Q\n5OG2z/mL/1AkSbMZZaX/WeAdU8Y2A1+vqt8A7gQuA0hyJvD6qjoOeC9wXRs/Cvgw8FbgZGDL8BOF\nJGlpzBr6VfVN4CdThs8FbmjbN7TLk+M3tv22A2uSrGXwpLGtqp6sqgPANmDjwsuXJM3FfHv6v15V\n+wGq6glgbRtfD+wdmvdYG5s6vq+NSZKW0OpFup2aYTzzubGJiYkXtnu9Hr1ebz43I0krVr/fp9/v\nz3m/VM2U10OTkqOBL1fVm9vl3UCvqvYnWQfcVVXHJ7mubX++zXsQ+C3gtDb/fW38F+ZNua8apSZ1\nTxJmXl+MdAv4s6WVKglVNetCe9T2TvjFVfutwIVt+0LglqHx81sBpwAHWhvoDuD0JGvaQd3T25gk\naQnN2t5JchPQA16R5AfAFuDjwH9N8i7gUeA8gKq6LclZSb4PPA1c1MZ/kuSjwLcYLNWuaAd0JUlL\naKT2zlKyvaOZ2N6RZrbY7R1J0gpg6EtShxj6ktQhhr4kdYihL0kdYuhLUocY+pLUIYa+JHWIoS9J\nHWLoS1KHGPqS1CGGviR1iKEvSR1i6EtShxj6ktQhhr4kdYihL0kdYuhLUocsKPSTPJLkO0nuS3JP\nGzsqybYkDyW5I8maofmfTLInyc4kJyy0eEnS3Cx0pf880KuqE6vqpDa2Gfh6Vf0GcCdwGUCSM4HX\nV9VxwHuB6xZ435KkOVpo6Gea2zgXuKFt39AuT47fCFBV24E1SdYu8P4lSXOw0NAv4I4k9yZ5Txtb\nW1X7AarqCWAy2NcDe4f23dfGJElLZPUC9z+1qh5P8mvAtiQPMXgiGDb1siRpTBYU+lX1ePv3b5P8\nFXASsD/J2qran2Qd8Ddt+j7g1UO7b2hjLzIxMfHCdq/Xo9frLaRMSVpx+v0+/X5/zvulan4L8SS/\nAqyqqqeSvATYBlwBvB34cVVdmWQz8E+qanOSs4CLq+rsJKcA11TVKdPcbs23Jq1sSVjYC8fgz5ZW\nqiRUVWabt5CV/lrgS0mq3c5/qaptSb4F3JzkXcCjwHkAVXVbkrOSfB94GrhoAfctSZqHea/0DxZX\n+pqJK31pZqOu9H1HriR1iKEvSR1i6EtShxj6ktQhhr4kdYihL0kdYuhLUocY+pLUIYa+JHWIoS8t\ngXXrjiHJvL/WrTtm3A9BK4R/hkGHjEP5zzAcyrXr0OCfYZBWlCN8paBF4Upfh4xDebW8GLUfqo9d\nS8OVviTpRQx9SeoQQ18a0ULOwJGWC0Nfc9LlUw/373+UQV99Pl/jNv8DwYfy/5lezAO5mpOFH5D8\nx8AzC9h/fAczF/bYx38gdyG1+zu5/C3FZ+RK8/AMCwtOSQux5O2dJBuTPJjk4SSXLvX9d91C2zOS\nDm1LGvpJVgF/AbwDeBPwu0l+cylrmK9+vz/uEl5kPjUtrC89ygp97jUdfP1xFzCN/rgLmEF/mrHx\nvjFspfzuLRdLvdI/CdhTVY9W1bPAVuDcJa5hXpbjf/JyrGl5hll/3AVMoz/uAmbQn2ZssiU3v6/B\nQmMBFS3Dn/PlWNOoljr01wN7hy4/1sY0B5MtmiuuuML2zJwsbMWq+VrY9/1P//SacT+AFeWQPGVz\nx44dC/ohSsKBAwfmff8L7YsfdthLFrT/z1s0Wzj0Th0cp4WtWDVfC/u+P/30U2P7XVuJp6su6Smb\nSU4BJqpqY7u8GaiqunJojr9dkjQPo5yyudShfxjwEPB24HHgHuB3q2r3khUhSR22pOfpV9VzST4A\nbGPQWrrewJekpbPs3pErSTp4lvWB3CR/lOT5JC9fBrV8JMl3ktyX5PYk68ZdE0CSq5LsTrIzyReT\nvGwZ1PRvkjyQ5LkkbxlzLcvqzYBJrk+yP8l3x13LpCQbktyZ5HtJ7k/yh8ugpiOSbG+/b/cn2TLu\nmiYlWZVkR5Jbx13LpCSPDOXTPb9s7rIN/SQbgNOBhZ3ku3iuqqp/VlUnAv+dwakzy8E24E1VdQKw\nB7hszPUA3A/8a+Ab4ywiy/PNgJ9t9SwnPwP+Q1W9CfgXwMXj/j5V1TPAae337QTgzCQnjbOmIZcA\nu8ZdxBTPA72qOrGqfun3admGPvBnwB+Pu4hJVfXU0MWXMPgmj11Vfb2qJmu5G9gwznoAquqhqtrD\n+P9YzrJ7M2BVfRP4yThrmKqqnqiqnW37KWA3y+D9M1X107Z5BIPjj2PvRbfF6FnAp8ddyxRhxDxf\nlqGf5Bxgb1XdP+5ahiX5WJIfAL8HfHjc9UzjXcBXx13EMuKbAecoyTEMVtbbx1vJC22U+4AngK9V\n1b3jromfL0bH/gQ0RQF3JLk3yb//ZRPH9lc2k3wNWDs8xKDwDwGXM2jtDF83zpo+WFVfrqoPAR9q\nveE/ACaWQ11tzgeBZ6vqpuVSkw4tSV4KfAG4ZMor27For2BPbMep/irJG6tqbG2VJGcD+6tqZ5Ie\n438lO+zUqno8ya8BX0uyu72qfJGxhX5VnT7deJJ/ChwDfCdJGLQrvp3kpKr6m3HUNI2bgNtYotCf\nra4kFzJ4yfm2pagH5vS9Gqd9wGuGLm9oY5oiyWoGgf+5qrpl3PUMq6q/S3IXsJHx9tJPBc5JchZw\nJPCrSW6sqvPHWBMAVfV4+/dvk3yJQWtz2tBfdu2dqnqgqtZV1euq6rUMXpKfeLADfzZJjh26+E4G\nfc+xS7KRwcvNc9rBr+VmnKuhe4Fjkxyd5HBgE7AczrgIy2uVCPAZYFdVXTvuQgCSvDLJmrZ9JINX\n/g+Os6aquryqXlNVr2Pws3Tncgj8JL/SXqWR5CXAGcADM81fdqE/jWJ5/IJ8PMl3k+wE/hWDI/jL\nwZ8DL2Xwkm5Hkk+Nu6Ak70yyFzgF+EqSsRxnqKrngMk3A34P2DruNwMmuQn438AbkvwgyUXjrKfV\ndCrw+8Db2il/O9piYpxeBdzVft+2A3dU1W1jrmm5Wgt8sx3/uBv4clVtm2myb86SpA45FFb6kqRF\nYuhLUocY+pLUIYa+JHWIoS9JHWLoS1KHGPqS1CGGviR1yP8H9yiXHBw7/LYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f52fbf99e50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(eta,20);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEACAYAAACwB81wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAFjlJREFUeJzt3X2wXVV98PHvLxAgLyQgkYBEiEQRI7TgC/QRqkcoEC1I\nnIFB+3SQOrbSKjqWUR5E6/UZR0tLsfL4Unmxw4sKlhJFBUlSPVARCEoCCCGACZCQhGCDIS8Q8rKe\nP/YOuYTchHvPOufse/b3M3Mn++6781trr8n53ZW1114rUkpIknrfiG5XQJLUGSZ8SaoJE74k1YQJ\nX5JqwoQvSTVhwpekmtg1R5CIeAxYBWwGNqSUjsoRV5KUT5aET5HoGymlZzLFkyRllmtIJzLGkiS1\nQa4knYBbIuLuiPjrTDElSRnlGtI5JqW0LCJeDcyKiPkppV9mii1JyiBLwk8pLSv/fDoiZgBHAS9J\n+BHhoj2SNAQppcgRp+UhnYgYHRFjy+MxwInAb7d3bUrJr5T4whe+0PU6VOXLtrAtbIsdf+WUo4c/\nEZhR9uB3Bb6bUpqZIa4kKaOWE35KaRFwRIa6SJLayKmUXdBoNLpdhcqwLbayLbayLdojco8RDVhQ\nROpUWZLUKyKCVJWHtpKk9njhhbzxTPiSVFFve1veeCZ8Saqo3KPgJnxJqgkTviRVlD18SaoJE74k\naUhM+JJUUfbwJakmTPiSpCEx4UtSRdnDlyQNiQlfkirKHr4k1URlE35EjIiIeyLixlwxJUn55Ozh\nfxJ4MGM8Saq1SvbwI2IS8F7g8hzxJEkVTfjAV4FPA25pJUkV1fIm5hHx58BTKaV5EdEABtyKq6+v\n78XjRqPhvpWStI1ms0mz2QTgf/4nb+yW97SNiC8DfwlsBEYBewI3pJTO3OY697SVpEGYMgUWLsy3\np23WTcwj4l3AuSml923nZyZ8SRqEgw+GRYvcxFySel7uPnLLY/j9pZRuBW7NGVOSlIc9fEmqqKpO\ny5QkZWbClyQNiQlfkirKHr4kaUhM+JJUURs25I1nwpekilq3Lm88E74kVVBKsHZt3pgmfEmqoJUr\nYY898sY04UtSBS1dCpMn541pwpekClqxAvbdN29ME74kVdCiRbDffnljmvAlqYIWLoSpU/PGNOFL\nUgWtXg3jx+eNacKXpAp69lnYc8+8MU34klRBzz4L48bljZljE/PdgduA3cp416eUvthqXEmqs9Wr\n8/fwW074KaX1EfHulNK6iNgFuD0ibk4pzclQP0mqpXYk/CxDOimlLSs+7E7xS8TdyiWpBWvXwpgx\neWNmSfgRMSIi5gLLgVkppbtzxJWkumpHws+yiXlKaTNwZESMA34YEVNTSg9ue11fX9+Lx41Gg0aj\nkaN4SeoZzWaTZrPJihXwrW/ljR0p85YqEfF5YG1K6eJtzqfcZUlSL1q0CA4+GJ57DkaNClJKkSNu\ny0M6ETEhIsaXx6OAE4CHWo0rSXV1223w9rfnXy0zx5DO/sCVETGC4hfIdSmlmzLElaRaWrMG3vrW\n/HFzTMu8H3hLhrpIkoBVq2CvvfLH9U1bSaqYBQvyr6MDJnxJqpx58/JvfgJtmKUzYEHO0pGkndq4\nEUaOhKefhgkTIKJCs3QkSfnMnFn8OWFC/tgmfEmqkMsug7/7u/bEzvKmrSQpj//6L7jllvbEdgxf\nkipi8+biZavVq2H33YtzjuFLUg/6xjdgw4atyT43E74kVcCmTfCJT8C//Vv7yjDhS1IF/PM/w9ix\n8NGPtq8ME74kddnq1TBnDlx0UXvLMeFLUpeddho88AAcfXR7y3GWjiR10eLFcOCBxQ5Xo0e//Oc5\nZ+mY8CWpi04/HZ55BmbP3v7PcyZ8X7ySpC5ZswbuuguuvLIz5TmGL0ldctJJsHIlHHlkZ8preUgn\nIiYBVwETgc3AZSmlS7ZznUM6klTatAl23RWeegr23Xfg6yo1hh8R+wH7pZTmRcRY4DfAqSmlh7a5\nzoQvSRSbkzcaMGIE3HHHjq+t1NIKKaXlKaV55fEaYD5wQKtxJalX/eQnMGrUwA9q2yXrQ9uImAwc\nAdyVM64k9YJVq4okf+GFcM45MGZMZ8vPlvDL4ZzrgU+WPf2X6evre/G40WjQaDRyFS9JlXf22fDo\no/C2t8EZZ2z/mmazSbPZbEv5WebhR8SuwE+Am1NKXxvgGsfwJdXW+vXFWjm//jX88R+/8r9XqYe2\nABFxFfD7lNLf7+AaE76k2tm4ES6/HG68sXjB6le/ghhE+q7UQ9uIOAb438BxETE3Iu6JiGmtV02S\nhr+PfQwuvhiOOAJ++MPBJfvcXFpBktokJZg4EWbMgGOOGVoMl1aQpApbtQoefrjYm3bUqOIhbRWY\n8CUpo40b4c/+rFj9ctw4+I//aN+WhYNlwpekjO69F5Ytg4ULYbfdul2bl3LxNEnK6JJL4Ljjqpfs\nwR6+JLVk1ixYtKg4Xr4cbrhh6/dVY8KXpCFICc49F665BqZP33r+mmtgwoTu1WtHnJYpSUPw4INw\n4olw660wZUr7yqnUi1eSVDfPPQdnnVVsT9jOZJ+bCV+SBunCC2HkSPjiF7tdk8FxDF+SBrBpU5Hc\nb7rppefvvRfmzCnm2Q8nJnxJ2o4HHoBTTim2IfzGN4o3ZrcYPRre9Kbu1W2ofGgrSdtx3nmwYUOx\n8Fk3uZaOJLXJypXFEM63v10sadxLTPiSVFq9Gk46CV71KviLv4B3vKPbNcrLhC9JFD37Y4+F17ym\n2GR85Mhu1yi/XDteXQGcDDyVUvqjAa5xDF9S17zwAsydC7fdVmwzuK0nn4R99oEf/ajzdduRKm5x\neCywBrjKhC+pik45pVjj5oADihem9tzz5dccf3z1lkWo3EPblNIvI+KgHLEkKbdHH4Vf/AJWrCim\nVNaVb9pK6mlLlxY7Tp19dr2TPWSch1/28H/skI6kKli6FL7+dbjvvmI+/U03wS67dLtWg1e5IZ1X\nqq+v78XjRqNBo9HoZPGSamLL9Mr994f3vrcYsx8uyb7ZbNJsNtsSO2cPfzJFD//wAX5uD19S2z35\nJEybBpMmwU9/CiOG+cB15ZZHjojvAb8CDomIJyLir3LElaTB+uxn4cgjeyPZ5+ZaOpJ6xpw5xdux\nTz8Ne+/d7drkMWzH8CUpp4UL4ctfhrVrYckSuPtu+Nd/7Z1kn5s9fEnD0saNcNppMH588YB21Cj4\n0z+t3otTrbKHL6n2vvSl4gHtd78LY8Z0uzbDgz18ScPOQw/B298Ov/kNHHJIt2vTXpWbpSNJnXLp\npcVuU+ef3/vJPjeHdCQNG1ddBf/wD/Db38Kb39zt2gw/DulIGhaeeAKmToVms1gbpy4qtzzyKyrI\nhC+pBe9/Pxx6KHzlK92uSWeZ8CXVypw5xXIJixfXb0aOD20l1cqMGfA3f1O/ZJ+bPXxJlbZkCRx2\nWPEW7Rve0O3adJ49fEm18e//DmecUc9kn5sJX1KlzZwJ06d3uxa9wSEdSZU1f37xRu3y5TB2bLdr\n0x0O6UiqhXvugZNPrm+yz82EL6my1qyBceO6XYvekWvHq2kR8VBEPBwR5+WIKUnPPw977NHtWvSO\nlhN+RIwAvg6cBLwZ+GBEHNpqXEl6/vlinXvlkaOHfxTwSErp8ZTSBuBa4NQMcSXV3HPP2cPPKUfC\nPwBY3O/7JeU5SWqJQzp59fTyyJ/5DFx9NUSLE5pyzCY1RnviGCN/jFxxcsR4/nn45jdbj6NCjoT/\nJHBgv+8nledepq+v78XjRqNBo9HIUPzArrwSbrwRXvva1mO1+kuj12LkimOM9sTppRh1m5LZbDZp\nNpttid3yi1cRsQuwADgeWAbMAT6YUpq/zXUdffFq/XrYc8+ihzDCyaeShqlKbWKeUtoUER8HZlI8\nE7hi22TfDcuXw8SJJntJ2iLLGH5K6WfAG3PEymXhQpg0qdu1kKTq6Nn+75IlMGVKt2shSdXRswn/\nd7+DAw/c+XWSVBc9OS1z5Ur49reLGTqSpEJP9vCvvhpOOKFYVlWSVOi59fCfew722guuvx5OOaXt\nxUlSW+WcltlzCf/jH4cVK+AHP2h7UZLUdib8ATzwABx9NDz2GEyY0NaiJKkj3PFqABdcAJ/6lMle\nkrZn2PfwFy+Gvj5Ytw6uvRYef9zpmJJ6h0M6/XzkI8U0zOnTYfJkeOc7sxchSV1TqbV0uiUluOgi\nuOIKeOKJPCtiSlIvG7Zj+P/yL3D55XD77SZ7SXolhm3Cv+EG+OpX4R3v6HZNJGl4GJYJ/2MfK6Zg\ntnn/FEnqKcPuoe2jj8KRR8KcOfCmN2WomCRVWG3n4a9cWbxYde65JntJGqyWEn5EnBYRv42ITRHx\nllyVGsjtt8Phhxfz7iVJg9NqD/9+4P3ArRnqslPXXgvTpnWiJEnqPS3Nw08pLQCIyLE3/fbiF0sd\n33EH3H9/8XXJJe0oSZJ6X6XH8JcuhQ9/GKZOhfPPh3nzYJ99ul0rSRqedtrDj4hZwMT+p4AEXJBS\n+vFgCuvrN/jeaDRo7GRe5UUXFVMvzzlnMKVI0vDVbDZpNpttiZ1lWmZE/AI4N6V0zw6uGdS0zEsv\nLVa//PWv4aCDWq6iJA1LVV1LJ9s4/n//d/Fy1e23m+wlKZdWp2VOj4jFwJ8AP4mIm1ut0JIl8PnP\nF8smHHVUq9EkSVtU6k3b55+H/feHY4+F666D0aM7UjVJqqyeXA9/xYpiTfuJE2HGjI5USZIqrycT\n/gc+AJs2wfe/D7sO21X6JSmvnkv4S5fClCmwfDmMH9+R6kjSsNBTi6etXg0nngjHH2+yl6R26nrC\n/9nPIAKuv77bNZGk3tb1hP+d78Bxx8Eee3S7JpLU27o6hj9/frFOzrPPwp57dqQakjSs9MxD2w9+\nENavL/anlSS9XFWXVhiUWbOK9e2XLOlWDSSpXro2ht9sFqtgHnBAt2ogSfXSlSGdlSuLde1nzy6m\nY0qStm/Yz8O/+WY45RSTvSR1UlcS/qxZ8J73dKNkSaqvriT8mTPhXe/qRsmSVF8dH8Nft65YQuGF\nF4o3bCVJAxvWY/hf+Qq8+tUme0nqtFZ3vPqniJgfEfMi4j8jYtzO/s7s2XDZZa2UKkkailZ7+DOB\nN6eUjgAeAc7f2V945BF461tbLFWSNGgtJfyU0uyU0uby2zuBSTu6/ve/h40bi12tJEmdlXMM/8PA\nDjcxX7AA3vhGx+8lqRt2upZORMwC+vfJA0jABSmlH5fXXABsSCl9b0exLr64j7Vroa8PGo0GjUZj\nyBWXpF7UbDZpNpttid3ytMyIOAv4a+C4lNL6HVyXLr00MWeOD20l6ZWqzGqZETEN+DTwzh0l+y1W\nrYIxY1opUZI0VK2O4f8/YCwwKyLuiYhv7ujim2+G0aNbLFGSNCQt9fBTSm8YzPVr1sDJJ7dSoiRp\nqDr6pu2yZbD//p0sUZK0RUfX0tljj8SqVbDbbh0pUpKGvWG7ls5rXmOyl6Ru6WjCf93rOlmaJKm/\njib8Qw7pZGmSpP46mvBf//pOliZJ6q+jCX/KlE6WJknqr6MJ/7DDOlmaJKm/jk7LfOGFxMiRHSlO\nknrCsJ2WabKXpO7p+J62kqTuMOFLUk2Y8CWpJkz4klQTJnxJqgkTviTVREsJPyL+b0TcGxFzI+Jn\nEbFfropJkvJq6cWriBibUlpTHp8DTE0p/e0A16ZOveQlSb2iMi9ebUn2pTHA5taqI0lql5b2tAWI\niC8BZwJ/AN7dco0kSW2x04QfEbOAif1PAQm4IKX045TS54DPRcR5wDlA30Cx+vq2/qjRaNBoNIZU\naUnqVc1mk2az2ZbY2RZPi4jXAjellA4f4OeO4UvSIFVmDD8i+m9pMh2Y31p1JEnt0uoY/j9GxCEU\nD2sfB85uvUqSpHbo6Hr4DulI0uBUZkhHkjR8mPAlqSZM+JJUEyZ8SaoJE74k1YQJX5JqwoQvSTVh\nwpekmjDhS1JNmPAlqSZM+JJUEyZ8SaoJE74k1YQJX5JqwoQvSTWRJeFHxLkRsTkiXpUjniQpv5YT\nfkRMAk6g2PFKr0C7NigejmyLrWyLrWyL9sjRw/8q8OkMcWrDf8xb2RZb2RZb2Rbt0eom5u8DFqeU\n7s9UH0lSm+x0E/OImAVM7H8KSMDngM9SDOf0/5kkqYKGvIl5RBwGzAbWUST6ScCTwFEppRXbud4d\nzCVpCHJtYj7khP+yQBGLgLeklJ7JElCSlFXOefgJh3QkqbKy9fAlSdXW9jdtI2JaRDwUEQ9HxHnt\nLq8bIuKKiHgqIu7rd27viJgZEQsi4paIGN/vZ5dExCMRMS8ijuh3/kNlOy2IiDM7fR85RMSkiPh5\nRDwQEfdHxCfK87Vrj4jYPSLuioi5ZVt8oTw/OSLuLO/t+xGxa3l+t4i4tmyLOyLiwH6xzi/Pz4+I\nE7t1T62KiBERcU9E3Fh+X8u2iIjHIuLe8t/GnPJc+z8jKaW2fVH8QnkUOAgYCcwDDm1nmd34Ao4F\njgDu63fuQuAz5fF5wD+Wx+8BfloeHw3cWR7vDfwOGA/steW42/c2hLbYDziiPB4LLAAOrXF7jC7/\n3AW4s7zH64DTy/PfAj5aHv8t8M3y+Azg2vJ4KjCXYlbd5PIzFd2+tyG2x6eAa4Aby+9r2RbAQmDv\nbc61/TPS7h7+UcAjKaXHU0obgGuBU9tcZsellH4JbPuw+lTgyvL4Srbe96nAVeXfuwsYHxETgZOA\nmSmlVSmlPwAzgWntrntuKaXlKaV55fEaYD7FDK66tse68nB3iiSVgHcD/1mevxKYXh73b6PrgePK\n4/dRJLyNKaXHgEcoPlvDSvlW/nuBy/udPo4atgXF885t82/bPyPtTvgHAIv7fb+kPFcH+6aUnoIi\nCbL1XYaB2mTb808yzNsqIiZT/M/nTmBiHdujHMKYCywHZlH0wv6QUtpcXtL/M/HiPaeUNgGryvWp\neqIt2PpWfgKIiH2AZ2raFgm4JSLujoiPlOfa/hnZ6YtXymagp+M9ObMpIsZS9Mw+mVJas533MGrR\nHmUyOzIixgEzKIa3XqmeaYuI+HPgqZTSvIho9P/RKw2Rv1ZddUxKaVlEvBqYGRELePlnIvtnpN09\n/CeBA/t9v+XlrDp4qvxvFxGxH7DlZbQngdf2u25Lm/RMW5UP3q4Hrk4p/ag8Xdv2AEgpPQs0gf8F\n7BURWz57/e/rxbaIiF2AcSmllQzcRsPJMcD7ImIh8H2KIZqvUQxP1K0tSCktK/98GvghxbBU2z8j\n7U74dwOvj4iDImI34APAjW0us1uCl/7mvRE4qzw+C/hRv/NnAkTEn1D89/4p4BbghIgYHxF7UyxZ\ncUv7q90W3wEeTCl9rd+52rVHREzYMtMiIkZR3MODwC+A08vLPsRL2+JD5fHpwM/7nf9AOXPldcDr\ngTntv4N8UkqfTSkdmFI6mCIP/Dyl9JfUsC0iYnT5P2AiYgxwInA/nfiMdOBp9DSKmRqPAP+n20/H\n23SP3wOWAuuBJ4C/oniCPru895nAXv2u/zrF7IJ7Kd5O3nL+rLKdHgbO7PZ9DbEtjgE2UczImgvc\nU/4beFXd2gM4vLz/ecB9wAXl+dcBd5X3dR0wsjy/O/CD8p7vBCb3i3V+2UbzgRO7fW8ttsu72DpL\np3ZtUd7zls/H/VvyYic+I754JUk14RaHklQTJnxJqgkTviTVhAlfkmrChC9JNWHCl6SaMOFLUk2Y\n8CWpJv4/xVJueuuIZhYAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f52f757ac10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(sorted(eta));"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
